Temperature dependent surface relaxations of Ag(lll) 
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The temperature dependent surface relaxation of Ag(lll) is calculated by density-functional theory. 
At a given temperature, the equilibrium geometry is determined by minimizing the Helmholtz 
free energy within the quasiharmonic approximation. To this end, phonon dispersions all over the 
Brillouin zone are determined from density-functional perturbation theory. We find that the top- 
layer relaxation of Ag(lll) changes from an inward contraction (—0.8%) to an outward expansion 
(+6.3%) as the temperature increases from T = K to 1150 K, in agreement with experimental 
findings. Also the calculated surface phonon dispersion curves at room temperature are in good 
agreement with helium scattering measurements. The mechanism driving this surface expansion is 
analyzed. 
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I. INTRODUCTION 

The equilibrium geometry of a system depends on the 
temperature due to the anharmonicity of the interatomic 
potential. The presence of a surface breaks the peri- 
odic structure normal to the surface, and anharmonic 
effects are expected to be larger at the surface than in 
the bulk [0. Hence, the surface interlayer separation may 
change more strongly with temperature than the bulk 
lattice parameter. Indeed, enhanced anharmonic effects 
have been observed by recent experiments on several sur- 
faces: NifOOl) §, Pb(llO) |, CufllO) §, Ag(lll) §, 
Cu(lll) @ as well as Be(0001) 0. Among them, the 
large thermal expansion observed in the close-packed 
Ag(lll) surface has attracted much attention P~p!o|, but 
at present the interpretation of these results is still con- 
troversial. 

Using an the embedded-atom method (EAM) in which 
the parameters of the interactomic potential are deter- 
mined by fitting bulk properties, Lewis j^] simulated the 
thermal behavior of Ag(lll) for a large range of temper- 
atures using molecular dynamics. The results for the 
top layer relaxation differ significantly from those re- 
ported by an experimental study Q: the top interlayer 
spacing, dx%, remains smaller than the bulk value even 
at temperatures as high as 1110 K; while the analysis 
of experimental results ||] obtained by the medium en- 
ergy ion scattering (MEIS), concluded that d\2 changes 
from —2.5% contraction to 10.0% expansion as tem- 
perature increases from room temperature to 1150 K. 
Narasimhan and Scheffler || investigated the tempera- 
ture dependence of d\2 by minimizing the Helmholtz free 
energy of the system with respect to d\2 in a simpli- 



fied quasiharmonic approximation (QHA), where the vi- 
brational free energy was calculated including only three 
representative modes corresponding to the rigid vibra- 
tion of the top layer on a rigid substrate. The static 
total energy and the vibrational frequencies, were calcu- 
lated using density-functional theory (DFT) within the 
local-density approximation (LDA) . Although the results 
obtained within this "three-mode approximation" over- 
estimated the effect (e.g., at T = 1040 K, the calcu- 
lated surface relaxation is 15% whereas the experimental 
value was 7.5% ), these calculations provided a physi- 
cal explanation of the mechanism underlying the thermal 
expansion observed at this surface. Subsequently, using 
again EAM potential, Kara et al. [jlO| obtained a rather 
small thermal expansion. They argued that the large 
thermal expansion of Ref. Q was the result of an im- 
proper representation of the vibrational density of states. 
On the other hand, very recent MEIS measurements on 
Cu(lll) § and LEED measurements on Be(0001) § 
seem to support the theoretical picture developed in Ref. 

!■ 

Recent calculations of the thermal properties of Ag 
bulk |ll| demonstrate that the QHA provides a very ac- 
curate description of the thermal expansion and heat ca- 
pacity of Ag up to the melting point. In order to resolve 
the controversy on the thermal behavior of Ag(lll), we 
have recalculated the surface thermal expansion of this 
surface within DFT-LDA and QHA without any further 
approximations. In particular, the vibrational contribu- 
tions to the free energy from the whole BZ are included 
thanks to the efficient calculation of phonon dispersions 
by density-functional perturbation theory fOg] . Our re- 
sults positively indicate that DFT and the QHA — at 
variance with previous attempts based on EAM [p|]ic|]— 
provide a quantitatively correct description of the anoma- 
lous thermal properties of this surface. Our results also 
show the importance of a proper sampling of vibra- 
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tional modes over the BZ for a quantitatively reliable 
result. The qualitative explanation of the earlier work of 
Narasimhan and Scheffler || is fully confirmed. The dis- 
agreement with reported EAM results [^|jic|l is argued to 
due to the approximate nature of EAM and to an incor- 
rect k-point summation in Ref. fl(i[|. 



II. COMPUTATIONAL DETAILS 

To model the surface, we adopt a repeated-slab ge- 
ometry consisting of seven atomic layers separated by a 
vacuum region corresponding to five atomic layers. As in 
a previous treatment p^ , the Helmholtz free energy of 
the slab is given by 



F({d},T) = E({d})+F vih ({d},T) 
= E({d}) + 

(hu) p {(\\u{d}) 



3N , 

k B T^2Yl lns2sinh 

qii p=l 



(2.1) 



where E is the static total energy, as obtained by DFT 
calculations, k B and h are the Boltzmann and the Planck 
constants, and {d} represents the set of interlayer dis- 
tances normal to the surface and the inter-particle dis- 
tances parallel to the surface. The vibrational free energy 
is denoted as F v ib, and w P (q|| , {d}) is the frequency of the 
p-th mode for a given wave vector q|| , evaluated at the ge- 
ometry defined by {d}; and N is the number of atoms in 
the slab. The static total energy E in Eq. (2.1) includes 
all the anharmonic terms of the interatomic potential. 
The anharmonic nature also appears in the vibrational 
free energy F v fo since in the quasiharmonic approxima- 
tion the vibrational frequencies ui p (q\\, {d}) are allowed 
to change with {d}. 

At a given temperature T and zero pressure, the 
equilibrium geometry is determined by the minimum 
of the Helmholtz free energy, i.e., dF/dd = dE/dd 
+dF V ih/dd — 0. It is not practical to minime the 
Helmholtz free energy with respect to all the lattice pa- 
rameters {d} within the present ab initio approach. We 
therefore assume here that the vibrational free energy 
F v ib depends parametrically only on the interlayer dis- 
tance between the top layers, d\2- All other interlayer 
distances and inter-particle distances are assumed to vary 
with temperature as in the bulk jHJ . The equation which 
determines the temperature dependent equilibrium di2, 
then reads 



dE{d 1 
ddi2 



dF v 



dd u 



dhuj p (q) 



(2.2) 



where n p is the occupation number of the p-th mode de- 
fined by 
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The static total energy E is calculated by density- 
functional theory within the local-density approxima- 
tion Hi . Fully separable norm-conserving pseudopoten- 
tials |jl5|] are used in our calculations together with a 
plane wave basis set with a kinetic energy cut-off of 55 
Ry. BZ integrations are performed with the smearing 
technique of Ref. |16| using the Hermite-Gauss smearing 
function of order one, a smearing width a = 50 mRy, 
and a 16-point grid in the irreducible wedge of the BZ. 
The phonon frequencies of the system are calculated by 
density- functional perturbation theory |l2] ] . The dynam- 
ical matrices are calculated on a 4x4 grid of points in 
the surface BZ of the 7-layer slab and Fourier interpo- 
lated over a 48-point grid of q| | vectors in the irreducible 
wedge of the surface BZ, in order to calculate the vibra- 
tional contribution to the free energy. 



III. RESULTS 

We first calculated the surface relaxation of Ag(lll) 
by minimizing the static total energy and neglecting the 
vibrational contributions. The obtained results are : 
Ad 12 /d B = -1.0%, Ad 23 /d B = -0.2% , where d B is the 
interlayer spacing in the bulk. Starting from this static 
equilibrium geometry, the Helmholtz free energy of the 
slab is evaluated as a function of the top layer interspac- 
ing d\2- The temperature dependence of the equilibrium 
d\2 is then obtained from Eq. (2.2). All other interlayer 
spacings are assumed to expand according to the tem- 
perature dependence of the bulk. 
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FIG. 1. Temperature dependence of surface layer relax- 
ation of Ag(lll). Our calculated results are denoted by the 
filled circles. Open circles are the results of Ref. ||, open tri- 
angles are EAM simulations Jl(| . The experimental results |^] 
are shown by filled diamonds. 

Figure 1 shows the calculated results (filled circles) to- 
gether with the experimental data S (filled diamonds) 
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and other calculations [p|jL0|]. The temperature is scaled 
with respect to the experimental melting temperature 
(T m = 1234 K ). It can be seen that in the low 
temperature range (T/T m < 0.6), the present calculation 
gives a thermal expansion similar to that of Ref. M . At 
higher temperatures, the calculations of Ref. || overes- 
timate the thermal expansion of Ag(lll). The present 
calculation, which includes all the vibrational modes of 
the slab in the whole BZ, displays a much smaller top 
layer expansion than obtained in the "three-mode ap- 
proximation" g], thus bringing the theoretical predic- 
tions in much closer agreement with experimental re- 
sults: Adi2/d,B changes from —0.8% inward contraction 
(including zero-point vibrations) at T = K to +6.3% 
outward expansion at T — 1150 K (the corresponding ex- 
perimental figures vary from —2.5% to +10.0% ||). The 
EAM simulation in Ref. flC|| , on the contrary, shows no 
enhancement of thermal expansion of the interlayer spac- 
ing in the whole temperature range, i.e., the di2 value 
always remains smaller than the interlayer spacing in the 
bulk. 




FIG. 2. Variation of the static total energy E, vibrational 
free energy _FVib and the Helmholtz free energy with the sur- 
face interlayer spacing d\2 at T=500 K. 

What is the mechanism giving Ag(lll) such an en- 
hanced thermal expansion? Figure 2 shows the variation 
of the static energy, E, the Helmholtz free energy, F, 
and the vibrational free energy, F v n,, with the top-layer 
interspacing di2, for T — 500 K. The equilibrium geom- 
etry, di2, is determined by two factors: one is the static 
total energy, E, which governs the binding strength of 
the surface with the substrate and the anharmonicity of 
the interlayer potential normal to the surface; the other 
is the decrease of the vibrational free energy, F v ib , which 
reflects the "softening" of the vibrational frequency with 
the increase of d\2- We note that does not always 



decrease with the increase of lattice parameters, for ex- 
ample, in the anomalous thermal expansion of bulk sil- 
icon pqj. As d\2 increases from the static equilibrium 
value, ^^2, the interlayer potential increases while the vi- 
brational free energy decreases, determining a new equi- 
librium spacing d\2, larger than cq 2 . 
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FIG. 3. Variation of —dF v -^/dd\2 (solid lines) and 
dE/dd\2 (dashed line) with d\2- The harmonic results of 
dE/dd\2 are shown by the dot-dashed line. The equilibrium 
di2 is obtained by the intersection of the solid lines and the 
dashed line. 

Figure 3 shows how the temperature dependent equi- 
librium d\2 is obtained from Eq. (2.2). The dashed line is 
the variation of dE/ddu with d\2- The dot-dashed line is 
the derivative of E with respect to d\2 when anharmonic 
terms in E are neglected by quadratically expanding it 
around the static equilibrium geometry. The solid lines 
are —dF vl b/ddi2 at different temperatures. The equilib- 
rium geometry of di2 is determined by the intersection of 
dEjddn and — dF^j dd\2- It can be clearly seen that by 
increasing the temperature, the intersection of dEjddyi 
and —dF v ib/ddi2 gives an increasing value of equilib- 
rium distance dyi- Furthermore, the anharmonicity of 
the static energy E, dashed line, is essential in determin- 
ing the enhanced thermal expansion at high temperature, 
that would be much reduced if this anharmonicity were 
neglected (dot-dashed line). 

The driving force for expansion comes from the tem- 
perature variation of the vibrational free energy. Equa- 
tion (2.2) reveals that the value of —dF^/dd\2 at a given 
temperature is determined by dhuj p /ddi2 which repre- 
sents the "softening" of oj p with the increase of d±2 ■ Such 
"softening" of the frequencies is characterized by the shift 
of the phonon density of states (DOS). 
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FIG. 4. The phonon density of states (DOS) of different 
atomic layers. Solid lines are for di2=2.32 A, dashed lines are 
for di2=2.42 A. 

Figure 4 shows the shift of layer localized phonon DOS 
with di2- The solid line corresponds to di2 = 2.32 A and 
the dashed line corresponds to d\2 = 2.42 A. Clearly, the 
DOS of the middle layer of the slab is quite bulk like, 
and the shift with d\2 is very small, while the surface- 
layer DOS is significantly different from that of the mid- 
dle layer and is very sensitive to d\2- By increasing d\2, 
the vibrational frequencies of the first layer "shift" down- 
ward, which results in a decrease of the vibrational free 
energy as shown in Fig. 2. The vibrations of atoms in 
the second layer also "soften" by increasing d\2 and con- 
tribute to the decrease of the vibrational free energy. 
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FIG. 5. The shift of phonon DOS at surface from differ- 
ent vibrational modes with different d\2- The solids are for 
di2=2.32 Aand dashed lines are for di2=2.42 A. 

Figure 5 shows the "shift" of the phonon DOS in the 
top layer with d\2 for different vibrational modes. It 
can be seen that as d\2 increases from 2.32 A to 2.42 
A, not only the vibrational modes that are perpendic- 
ular to surface, but also those parallel to surface, shift 
downward. The shift of the perpendicular modes reflects 
the anharmonicity of the interlayer potential E normal 
to the surface, while the shift of the parallel modes is 
due to the flattening of the corrugation of the poten- 
tial parallel to the surface. The latter effect implies that 
anharmonicity of the interatomic potential at a surface 
is a somewhat complicated effect where the x, y, and z 
degrees of freedom are not independent. The modes per- 
pendicular to the surface (corresponding to the Rayleigh 
wave) are mainly located in the low frequency range and 
provide a smaller contribution to the DOS. The parallel 
modes have relatively higher frequencies and occupy a 
larger portion of the total DOS. 

As discussed above in Fig. 2, the increase of 
—dF v ib/ddi2 with temperature determines the expansion 
of di 2 . 
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FIG. 6. Contributions to the vibrational free energy from 
different vibrational modes 

In Figure 6 the different contributions to —dFy^/ddvi 
coming from perpendicular and parallel modes are anal- 
ysed. It can be seen that the contribution of the parallel 
modes to —dF^/dd^ is larger than the one of perpen- 
dicular modes and that the difference between the contri- 
butions from parallel and perpendicular modes increases 
with temperature. This is in agreement with previous 
findings |9p5}]. Note, however, that this effect is not gen- 
eral and in other systems, as for instance the Be(0001) 
surface jl9| , the "softening" of parallel modes plays a mi- 
nor role. Altogether, the enhanced thermal expansion of 
di2 at Ag(lll) surface is governed by the anharmonicity 
of the interlayer potential normal to the surface as well 
as the "softening" of the parallel modes with the increase 
of di 2 . 




Finally we show in Figure 7 the surface phonon band 
structure corresponding to the geometry obtained for 
T = 300 K. The thickness of the slab has been extended 
to 30 atomic layers in order to decouple the surface vi- 
brations that penetrate deeply in the bulk, by inserting 
in the middle of our slab a number of layers with bulk- 
like force constants. The lowest surface-mode branch Si 
which lies below the bulk bands is the Rayleigh wave. 
It is primarily associated with vibrations normal to the 
surface ( shear- vertical (SV) mode, compare also Fig. 5). 
The gap mode S3 is primarily a "shear-horizontal" (SH) 
mode which is associated with vibrations in the direction 
transverse to q|| and parallel to the surface. The surface 
modes S2 and S4 are primarily longitudinal modes which 
are associated with vibrations along the direction parallel 
to q||. The calculated frequencies of the Rayleigh mode 
Si are in good agreement with experimental data from 
helium scattering [2C| ]. 

IV. CONCLUSION 

Our calculation of the thermal properties of Ag(lll) 
surface are in quantitative agreement with the enhanced 
thermal expansion found in the experiments. This be- 
havior is found to be determined by two effects: i) the 
anharmonicity of the interlayer potential normal to the 
surface and ii) the decrease of the vibrational free en- 
ergy with increasing dyi- The latter effect reflects not 
only the anharmonicity of the interlayer potential nor- 
mal to the surface, but also the flattening of the corru- 
gation of the interlayer potential parallel to the surface 
with the increase of dyi- A recent calculation using the 
EAM shows that the interlayer spacing dii increases 
very weakly between T = K and 1100 K. This sig- 
nificant difference to our ab initio result is due to two 
reasons: i) The k-summation in Ref. JTo[ ] was done in- 
correctly, i.e., the ratio of weights at f, K and M should 
be 1:2:3, not 1:6:6, and at the fee (111) surface the vi- 
brational modes at M and K are different, ii) the EAM 
lacks some important aspects of the quantum mechanical 
nature of electrons |21f| . 
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FIG. 7. Calculated phonon dispersion of a 30-layer slab 
modeling an Ag (111) surface. 
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